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1.        INTRODUCTION 

The  exhaust  of  a  large  space-based  HF/DF  chemical  laser  can  be  idealized  as  a  zero-thrust 
supersonic  ring- symmetric  jet  (Figure  1-1).  Assuming  a  vacuum  background,  the  exhaust  plume  is 
always  flanked  by  a  pair  of  ring-symmetric  rarefaction  waves  centered  at  the  nozzle  lips  (Figure  1-1). 
When  ambient  molecules  traveling  at  orbital  speed  impinge  obliquely  at  the  centered  rarefaction  wave 
(CRW),  they  give  rise  to  a  molecular  backflow  of  scattered  exhaust  species.  This  effect  constitutes  a 
potentially  significant  contribution  to  spacecraft  contamination  [1]. 

iMost  ambient  molecules  are  stopped  within  several  mean  free  paths  from  their  point  of  entry  into 
the  plume.  A  quantitative  estimate  of  ambient  back-scattering  would  thus  depend  on  the  flow  field  at 
the  outer  (hypersonic)  fringes  of  the  lip-centered  CRW.  Even  though  the  flow  in  those  regions  is 
generally  past  the  point  of  continuum  breakdown  [2],  the  density  there  is  reasonably  well 
approximated  by  the  continuum  flow  field,  as  demonstrated  by  Bird's  Monte-Carlo  simulation  of  a 
Prandtl-Meyer  expansion  to  vacuum  [3].  The  evaluation  of  ambient  scattering  thus  calls  for  an 
ancillary  computational  procedure  capable  of  rendering  the  continuum  flow  field  at  a  large  number  of 
points  in  the  ring-symmetric  CRW  of  an  HF/DF  laser  exhaust  plume. 

The  purpose  of  this  report  is  to  present  a  doubleheader  approach  to  this  CFD  task,  consisting  of  a 
specially  formulated  finite  difference  scheme  valid  throughout  the  plume  and  an  analytic 
approximation  for  the  CRW  portion  of  the  flow.  This  approach  is  motivated  by  [he  need  to 
approximate  the  CRW  tlovv  field  in  a  simpie  and  computationally  affordable  way.  A  finite  difference 
integration  by  marching  out  from  the  nozzle  exit  to  every  point  in  the  CRW  where  ambient  scattering 
is  to  be  evaluated,  is  not  affordable  due  to  the  very  large  number  of  such  points.  An  interpolation 
from  a  2-D  grid  of  pre-computed  points  is  affordable  but  rather  cumbersome  and  complex  (if  only  for 
the  need  to  maintain  an  elaborate  2-D  finite  difference  code  that  would  make  its  output  available  to 
the  ambient  scattering  code).  We  propose  to  obviate  both  the  finite  diiference  computation  and 
subsequent  interpolations  by  constructing  an  analytic  approximation  to  a  ring-symmetric  CRW.  The 
finite  difference  code  "JET"  that  was  written  for  the  purpose  of  computing  a  ring-symmetric 
supersonic  flow  field,  is  thereby  relegated  to  the  role  of  aiding  in  the  verification  of  the  analytic 
approximation.  We  aiso  use  this  code  to  obtain  whole-plume  solutions  for  the  purpose  of  illustrating 
some  features  of  the  flow  field  by  means  of  graphic  output. 

The  finite  difference  scheme  used  for  computing  the  exhaust  flow  field,  is  a  modification  of  the 
well-known    inverse    marching    characteristic    method    [4].      Rather    than    using    the    two    velocity 


components  often  recommended  in  the  literature  [4],  we  use  the  two  Riemann  invariants  (v±0)  as 
flow  variables.  The  key  element  in  the  scheme  is  a  Semi-Inverse  Marching  Algorithm  (SI MA).  The 
flow  is  assumed  to  exit  the  spacecraft  with  a  uniform  supersonic  speed  in  the  radial  (y)  direction, 
which  is  hence  the  marching  direction.  New  grid  points  are  determined  by  the  forward  intersection  of 
continuous  C+  (fan)  characteristics  with  a  new  line  y=ynew  .  whereas  segments  of  the  transverse 
C~  characteristic  lines  are  reversely  extended  from  each  new  grid  point  and  require  interpolation 
between  old  line  grid  points  (y  =  yoW).  This  modified  marching  scheme  is  analogous  to  existing  time 
dependent  1-D  characteristics  methods  (see  Sections  19-6(a)  and  19-6(j)  in  [4]  );  however,  we  found  no 
reference  to  the  use  of  Riemann  invariants  as  flow  variables  in  a  CRW  computation  in  order  to 
reduce  interpolation  errors. 

The  resulting  SIMA  scheme  replicates  a  planar  CRW  (Prandtl-Meyer  flow)  exactly,  lending  an 
extra  measure  of  credibility  to  its  accuracy  in  computing  a  ring-symmetric  CRW.  For  the  new  line 
segment  lying  outside  the  CRW,  grid  points  are  evenly  distributed  through  the  segment,  and  a  fully 
inverse  marching  scheme  is  used.  The  SIMA  scheme  is  described  briefly  in  chapter  2  ;  it  was 
implemented  in  a  FORTRAN  code  named  "JET",  which  was  specifically  written  for  the  computation 
of  ring-jets  having  a  vacuum  background.  A  detailed  description  of  the  code  JET  and  the  finite 
difference  schemes  on  which  it  is  founded  will  be  given  in  a  future  report  [1 1]  . 

Our  analytic  approximation '  to  a  ring-symmetric  CRW  is  formulated  as  follows.  In  a  planar 
(Prandtl-Meyer)  CRW  the  flow  is  uniform  along  the  characteristic  lines  that  fan  out  from  the  corner 
(we  assume  it  is  the  C  family).  In  the  ring- symmetric  case  the  flow  near  the  corner  approaches 
asymptotically  a  corresponding  planar  CRW  flow,  which  we  term  the  associate  CRW.  However, 
the  gradients  along  C  characteristics  at  the  corner  of  a  ring-symmetric  CRW  do  not  vanish  as  in  a 
planar  CRW.  The  key  idea  is  thus  :  evaluate  flow  gradients  in  C  "  directions  at  the  corner,  then 
use  them  to  extrapolate  the  associate  CRW  along  C  ^  lines  to  a  finite  distance  from  the  corner. 
This  extrapolation  constitutes  an  approximation  to  the  ring- symmetric  CRW.  Our  present  approach 
is  analogous  to  the  GRP  (Generalized  Riemann  Problem)  analysis  from  which  high  resolution  upwind 
schemes  for  time  dependent  Euler  equation  were  derived  [5,0]  .  In  terms  of  specific  results,  however, 
it  is  quite  different  from  the  original  GRP  analysis. 

A  perturbation  approach  to  steady  supersonic  flow  in  two-dimensional  space  (plane  or  axial 
symmetry)  using  characteristic  coordinates,  has  been  developed  by  Ostwatitsch  and  colleagues  [7]. 
Their  work  also  included  treatment  of  axisymmetric  CRW.  However,  their  approach  is  formulated  in 
terms  of  small  perturbation  relative  to  a  uniform  supersonic  flow.     In  our  analysis  the  CRW  is 


assumed  to  span  the  range  from  some  finite  pressure  to  vacuum,  and  the  perturbation  scheme 
consists  of  regarding  the  axisymmetric  terms  in  the  governing  equations  as  causing  a  small  deviation 
from  a  Prandtl-Meyer  flow.  Our  analytic  approximation  to  a  ring-symmetric  CRW  is  presented  in 
chapter  3  and  Appendix  A. 

The  present  approach  to  the  approximation  of  a  ring-symmetric  CRW  can  also  be  adapted  to  treat 
other  axisymmetric  centered  waves,  such  as  the  divergent  lip-centered  CRW  at  an  axisymmetric 
nozzle  exit,  or  the  cylindrically  converging  CRW  at  the  base  corner  of  an  axisymmetric  projectile 
moving  at  supersonic  speed  in  air. 

A  series  of  computations  were  performed  on  a  sample  case  of  typical  HF/DF  laser  exhaust.  A 
comparison  was  made  between  results  of  numerical  integration  (SI MA)  and  the  method  of  matched 
approximation  to  the  ring- symmetric  CRW  .  In  particular,  molecular  opacity  (expected  number  of 
collisions  along  a  path  of  a  fast  penetrating  molecule)  was  evaluated.  Reasonable  agreement  between 
SIMA  and  approximated  opacity  was  demonstrated.  These  results  and  their  analysis  are  presented  in 
chapter  4. 

It  seems  that  the  present  approach  can  be  adapted  to  other  problems  calling  for  opacity  of  the 
CRW  region  in  a  ring-symmetric  or  axisymmetric  exhaust  plume. 
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Figure    1-1. 


Ring-Symmetric   HFDF    Laser  Exhaust  Plume 


2.       THE   FINITE  DIFFERENCE  SCHEMES   FOR   COMPUTING   RING  -  JETS 

The  free  expansion  of  a  spacecraft  exhaust  plume  is  idealized  in  our  study  as  a  steady  supersonic 
isentropic  flow  in  two-dimensional  axisymmetric  coordinates,  with  a  zero-pressure  background.  The 
most  accurate  finite  difference  scheme  for  this  type  of  flow  is  the  well  known  method  of 
characteristics  [8].  At  an  earlier  phase  of  the  present  laser  exhaust  study,  a  code  AXSYM  [9]  was 
written  for  the  computation  of  ring-jet  flow  fields,  using  the  direct  method  of  characteristics.  A 
notorious  shortcoming  of  this  method  is  that  it  yields  the  flow  field  at  a  set  of  grid  points  formed  by 
the  (oblique)  intersection  of  the  C  and  C  families  of  characteristic  lines.  The  resulting  grid  is 
highly  irregular  (especially  at  regions  of  hypersonic  flow),  and  it  also  requires  the  retention  of  grid- 
image  matrices  for  the  several  regions  formed  by  intersection  of  families  of  characteristic  lines  and 
their  reflection  from  the  mid-plane  of  symmetry  (x  =  0). 

A  commonly  accepted  remedy  to  these  shortcomings  is  the  inverse  marching  characteristic  scheme 
[4].  The  marching  is  in  the  downstream  direction,  i.e.,  the  y  direction  in  our  case.  The  grid  points 
are  located  on  a  succession  of  constant  y  rows,  thereby  introducing  a  measure  of  regularity  in  the 
solution  grid.  Computer  memory  requirements  are  drastically  reduced  :  just  two  rows  are  kept  in 
core  memory  -  an  "old"  line  of  grid  points  and  a  "new"  line  of  grid  points.  For  reasons  which  will  be 
elucidated  below,  the  flow  variables  in  our  scheme  are  the  Riemann  invariants  ( v  ±  0)  .  The 
integration  of  of  the  How  equations  in  characteristic  form  (sometimes  referred  to  as  compatibility 
relations  [-4]),  is  performed  by  a  combination  of  two  marching  schemes.  At  grid  points  outside  the 
CRW  we  use  the  conventional  inverse  marching  scheme.  At  grid  points  within  the  CRW.  we  use  a 
modified  scheme  named  SI  MA  -  Semi  Inverse  Marching  Algorithm,  tailored  specifically  to  render 
accurate  computation  of  a  centered  wave  flow.  In  the  sequel,  we  outline  both  schemes  and  describe 
the  procedure  by  which  they  are  combined  to  yield  the  flow  field  of  a  ring-symmetric  jet.  More 
information  on  the  schemes  and  the  code  JET  will  be  provided  in  a  future  report  [1  lj  . 

The  basic  building  block  of  both  inverse  and  semi-inverse  marching  schemes  is  the  evaluation  of 
flow  variables  at  a  new  grid  point  x4  on  the  new  line  ty=ynew).  given  the  flow  at  a  row  of  grid 
points  on  the  old  line  (y=y  id),  which  is  initially  the  nozzle  exit  surface  where  flow  is  assumed 
uniform  (Figure  2-1). 

Consider  first  the  inverse  marching  scheme  (Figure  2-la).  The  trace  points  x;  .  x,  are 
determined  by  reversely  extending  C  and  C  '  characteristic  lines  from  point  (^4.ynew)  t0  the  old 
line  (Figure  2-la).    The  characteristic  segments  are  approximated  by  straight  line  segments,  whose 


slopes  are  initially  taken  from  linearly-interpolated  flow  properties  at  the  old  point  (x4,yold).  Trace 
points  X.  and  X,  are  assigned  values  of  flow  variables  obtained  by  linear  interpolation  between  their 
respective  nearest-neighbor  (old)  grid  points. 

The  compatibility  relations  along  segments  of  C  ~  and  C  characteristic  lines  in  finite  difference 
form  are  now  solved,  yielding  the  flow  variables  (v  ±  0)   at  the  new  grid  point  : 

Along  C+     (v-9)4  =  (v-8)2  +  co  sinn24  sinG24  An  /  y24 

(2-1) 
Along  C~     (v  +  9)4  =  (v  +  0),  +  co  sinfi,4  sin8,4  A^  /  y14 

Where  A^  ,  At]  are  the  length  of  the  respective  characteristic  segments;  indices  14,  24  refer  to 
centered  segment  values  obtained  by  averaging  the  values  of  variables  at  segment  endpoints.  The 
symmetry  index  co  is  as  follows :  co  =  0  for  plane  flow,  co  =  1  for  axisymmetric  flow.  The  usual 
isentropic  relation  [S]  is  used  to  determine  fi   from  v. 

Equation  (2-1)  is  now  regarded  as  an  implicit  relation  between  the  flow  variables  at  the  new  grid 
point  x4  and  the  interpolated  flow  variables  at  the  trace  points  X.  and  x.,  .  An  updated  pair  of 
trace  points  is  re-computed  from  an  updated  value  of  tlow  variables  at  the  new  grid  point,  and  the 
procedure  is  repeated  until  convergence  is  established. 

The  Semi  Inverse  Marching  Algorithm  (SIMA)  is  a  relatively  simple  modification  of  the  inverse 
marching  scheme.  Rather  than  seek  a  solution  on  a  new  grid  point  x4  whose  location  is  unrelated 
to  the  row  of  old  grid  points,  we  determine  x,  by  the  forward  extension  of  a  C'  characteristic  line 
from  an  old  grid  point  x.,  (Figure  2- lb).  The  trace  point  X,  is  determined  by  reversely  extending 
the  C  characteristic  line  from  the  new  grid  point,  just  as  in  the  inverse  marching  scheme.  The 
same  compatibility  equation  (2-1)  is  solved  for  the  flow  at  the  new  grid  point,  except  for  an  obvious 
geometrical  modification  :  whereas  in  the  inverse  marching  scheme  the  trace  points  X.  and  X,  were 
re-computed  in  each  iteration,  the  SIMA  variant  calls  for  re-computing  the  new  grid  point  x4  and 
the  trace  point   x{    until  convergence  is  established. 

The  resulting  scheme  replicates  a  (planar)  Prandtl- Meyer  flow  exactly.  The  reason  for  that  is  the 
combination  of  the  semi-inverse  marching  idea  with  the  choice  of  Riemann  invariants  as  flow 
variables.  Due  to  this  choice,  the  compatibility  relation  along  the  fan  characteristics  (C  )  reduces  to 
the  exact  relation    (v  -0)=  constant.    The  equation  along    C~      still  requires  interpolation  in  old 


values  of  (v  +  6).  However,  in  a  C  Prandtl- Meyer  flow  (v  +  0)  is  uniformly  constant,  so  that 
interpolation  does  not  introduce  any  truncation  errors.  It  is  noted  that  this  feature  is  lost  if  any  other 
flow  variable  is  used  in  conjunction  with  a  SI  MA  scheme  for  computing  the  flow  in  a  CRW. 

The  boundary  conditions  are  quite  simple.  Marching  starts  out  from  the  nozzle  exit  surface  where 
the  flow  is  assumed  uniform.  The  flow  is  bounded  on  the  left  by  a  mid-plane  of  symmetry  (x=0), 
where  the  boundary  condition  is  simply  9  =  90°.  On  the  vacuum  side  we  approximate  the  idealized 
zero-pressure  background  by  terminating  the  computation  at  a  high  Mach  number  (fan)  characteristic 
line  (typically  M  =  34  at  the  corner).  It  is  noted  that  as  a  result  of  ending  the  computation  at  a 
characteristic  line,  the  total  mass  flow  through  a  solution  line    y         decreases  slishtlv  as    v 

c  J  new  °      *  «  new 

increases. 

A  sample  computation  performed  by  the  code  JET  is  displayed  in  Figure  2-2.  The  CRW  region  is 
clearly  shown  as  bounded  by  the  final  C  characteristic  line  on  the  vacuum  side,  and  by  the  leading 
C  characteristic  and  its  reflection,  on  the  other  side.  The  code  JET  can  also  plot  an  assortment  of 
special  lines  :  characteristics,  continuum  breakdown  lines  (2],  lines  of  constant  Mach  number, 
streamlines  and  lines  of  constant  lateral  molecular  opacity  (expected  number  of  collisions  by  a  fast 
molecule  entering  the  plume  in  the  x  direction). 
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Figure   2-1.  Characteristic  Marching  Scheme  for  Supersonic  Flow  in  Two  Space  Dimensions 
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Figure   2-2. 


Rine-Svmmetric  Flow  Field  with  Centered  Rarefaction  Wave 


3.        ANALYTIC   APPROXIMATION   OF   RING-SYMMETRIC 
CENTERED   RAREFACTION  WAVES 

The  procedure  for  the  analytic  approximation  of  a  ring-symmetric  CRW  comprises  of  two  distinct 
phases.  The  first  phase  is  an  analysis  of  the  flow  field  near  the  corner  of  a  ring-symmetric  CRW, 
resulting  in  fan-wise  gradients  of  flow  variables  at  the  corner.  The  second  phase  is  an  extrapolation 
of  the  associate  CRW  at  the  corner,  using  the  fan-wise  gradients  obtained  in  the  analytic  phase.  A 
concise  description  of  each  phase  is  given  below,  followed  by  a  detailed  presentation  in  sections  3.1, 
3.2  and  3.3. 

The  first  phase  is  an  analysis  of  the  flow  near  the  corner  of  a  ring- symmetric  CRW.  It  is  noted 
that  in  a  planar  C  CRW  all  directional  derivatives  of  flow  variables  along  C  vanish,  but 
directional  derivatives  in  any  non  C  direction  increase  beyond  bound  as  the  point  at  which  they 
are  evaluated  approaches  the  corner.  This  suggests  that  the  analysis  of  a  ring-symmetric  CRW  as  a 
perturbation  of  its  associate  CRW  should  be  formulated  in  characteristic  coordinates.  The  goal  of 
this  analysis  would  be  to  derive  closed-form  expressions  for  the  C  directional  derivatives  of  flow 
variables  at  the  corner.   A  detailed  presentation  of  the  analytic  phase  is  given  in  section  3.1  below. 

In  the  second  phase,  the  C  ^  gradients  at  the  corner  are  used  to  extrapolate  the  associate  CRW 
at  the  corner  to  some  finite  distance  along  C  characteristic  lines.  It  is  this  extrapolation  that 
constitutes  our  analytic  approximation  to  a  ring- symmetric  CRW.  Rather  than  merely  extend  the 
associate  CRW  through  linear  extrapolation  (i.e..  a  Taylor  series  truncated  after  the  first-order  term) 
using  the  corner  C  derivatives,  we  opt  for  a  specially  formulated  "matched  extrapolation" 
scheme,  where  the  extended  tlow  field  is  matched  to  conform  exactly  to  the  flow  along  the  pre- 
reflection  segment  of  the  leading  characteristic  C"1"^)  .  The  result  is  a  ring-symmetric  CRW 
approximation  that  maintains  reasonable  accuracy  levels  up  to  several  corner  radii  away  from  the 
corner  -  a  considerably  larger  range  of  validity  than  that  of  a  linear  (unmatched)  extrapolation.  The 
matched  extrapolation  is  described  in  section  3.2  below. 

Finally  we  consider  an  "inverse  problem"  which  is  stated  as  follows.  Given  a  point  (x0,y0)  within  a 
ring- symmetric  CRW,  find  the  flow  variables  at  this  point.  The  matched  extrapolation  scheme  is 
geared  to  deliver  the  flow  field  along  entire  C+  characteristic  lines  (including  the  determination  of 
the  characteristic  lines  themselves  in  (x.y)  coordinates).  We  do  not  know  apriori  the  line  C  + 
passing  through  (xQ,y0)  .  and  we  seek  an  algorithm  for  obtaining  the  flow  at  (x0,y0)  without 
resorting  to  excessive  numerical  integrations.  The  solution  to  the  inverse  problem  relies  on  the 
observation  that   C      lines  are  usually  just  slightly  curved:  it  is  presented  in  section  3.3  below. 
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3.1      Analysis  of  Ring-Symmetric  CRW 

Consider  a  2-D  axisymmetric  steady  inviscid  flow  of  an  ideal  gas.  If  the  flow  is  also  supersonic, 
isentropic  and  irrotational,  the  governing  equations  in  characteristic  form  using  the  Riemann 
invariants  as  flow  variables  [8]   are  : 

Along  C  +  (p)    ....     (v  -  8)a  =  sinn(a,P)  sin9(a,P)  na(a,P)  /  y(a,P) 

Along  C  "  (a)    ....     (v  +  G)p  =  sinn(a,p)  sin6(a,p)  ^(a,p)  /  y(a,P) 

(3.1-1) 

Direction  of  C  +  (P)    V(a,P)  =  6(a,p)  +  n(a,P) 

Direction  of  C"(a)    X(«,P)  =  6(a,P)  -  n(a,P) 

where  (a,P)  are  the  characteristic  coordinates  ,  a  being  constant  along  C  ~~  (a)  and  P  being 
constant  along  C  (P).  ^(a.P)  and  r|(a,p)  are  the  distance  (from  an  as  vet  unspecified  origin) 
along    C  "     and    C    '  respectively.    We  note  that   (s,H)    would  not  qualify  as  an  alternate  set  of 

characteristic  coordinates  .    since  in  general    ^    is  not  constant  along    C  ^    and  likewise    x\    is  not 
constant  along   C     . 

The  specific  definition  of  (a.p)  ,  s(a-P)  an^  n(a-P)  1S  now  tailored  to  the  needs  of  the  intended 
ring-symmetric  CRW  analysis.  To  fix  ideas  we  consider  a  C  CRW  ;  it  is  then  natural  to  stipulate 
that  a  =  0  and  r|=0  at  the  corner.  Denoting  by  index  1  the  leading  characteristic  CT(p1j  .  we 
define  the  characteristic  coordinates   a    and   p   as  follows  : 

a  =  rua-P,) 

(3.1-2) 
P  =  M(O.P) 

Where  M(0,p)  is  the  Mach  number  at  the  associate  CRW.  We  shall  generally  use  p  in  lieu  of 
M(O.P)  ,  in  order  to  simplify  notation  while  maintaining  the  important  distinction  between  the 
associate  Mach  number   M(O.P)    and   the  Mach  number  within  the  ring-symmetric  CRW   M(a.p)  . 
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Let  Q(a,P)  denote  any  flow  variable.  We  are  seeking  the  corner  gradient  Qa(0,P)  in  a  ring- 
symmetric  CRW.  Since  any  flow  variable  Q  can  be  expressed  as  function  of  (v±0)  ,  Qa  can 
likewise  be  expressed  in  terms  of  (v±0)  and  (v±8)a  .  Consequently,  we  should  seek  to  derive 
expressions  for  (v±9)a  at  a  =  0  ,  from  the  governing  equations  (3.1-1).  These  equations  are  in  a 
particularly  simple  form,  in  that  only  one  Riemann  invariant  appears  in  each  equation.  Thus,  the 
corner  gradient  (v_6)a  is  already  given  by  (3.1-1)  upon  setting  a  =  0  .  We  now  turn  to  the 
derivation  of  (v  +  6)a  • 

A  closed  form  expression  for  (v  +  0)a  can  be  derived  from  (3.1-1)  as  follows.  Differentiate  the 
equation  along  C  characteristics  with  respect  to  a  ,  set  a  =  0  and  then  integrate  the  ensuing 
relation  with  respect  to  p  ,  from  the  leading  characteristic  C  (pj)  to  some  internal  characteristic 
C  (P)  .  Since  C  — (a)  characteristic  lines  within  the  CRW  shrink  in  size  as  a  approaches  zero, 
^p(0,p)  =  0  .  Hence,  the  only  non  zero  term  left  for  the  P  integration  along  C-(0)  is  the  term 
containing  ^ap(0,P)  .     The  result  is  the  following  expressions  for  (v  ±  G)a  at  the  corner  : 

[(v-e)J(0,P)  =  sinn(0,P)  sin9(0,p)  na(0.P)  /  y(0,P) 

(3.1-3) 


i-i 


P 
[(v  +  9)J(0.P)  =  shHKO.Pj)  /v(O.Pj)  +  [v(0.p,)j  '      j    sinu(O.m)  sin0(O.m)  s»(0.m)  dm 

Pi 


The  boundary  condition  for  the  integration  of  ( v  +  0)«    was  obtained  by  noting  that  along  the 
leading  characteristic  C     (p.)  : 

e(a.p,)  =  k/2 

ea(a.p,)  =  0  (3.1-4) 

(v-e)a  =  shnua.Pj)  /  y(a.P;) 


We  also  note  that  at  the  corner  V(O.p)  =  yc  for  any  p  .  However,  we  shall  maintain  the 
notation  y(a.p)  and  V(O.p)  in  order  to  emphasize  that  y(a,P)  is  a  field  variable.  The  geometrical 
derivatives  ^(O.p)  and  n.a(0.P)  are  readily  derived  from  standard  expressions  of  the  associate 
CRW  (Prandtl- Meyer  flow,  see  Appendix  A  for  details).   They  are  given  by  : 
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n„(o,P)  =  {[i+aY-iWl/Ii+aY-D/^,2]}^7'1"  {(Pl2- i)/(p2- 1)},/4 

(3.1-5) 


3  ,_2         x-1  r  2,-1 


Sap(0,P)  =  ((Y+D/4)no(0,P)P    (P   -l)"    [l+((Y-D/2)p] 


Now,  equations  (3.1-3)  and  (3.1-5)  constitute  closed  form  expressions  for  the  corner  gradients 
(v±6)a  .  The  only  additional  information  needed  is  the  standard  Prandtl-Meyer  expressions  for 
0(0, p)  and  v(0.p)  in  the  associate  CRW  (Appendix  A).  The  integral  in  (3.1-3)  cannot  be  expressed 
in  terms  of  elementary*  functions  ;  however,  it  is  readily  evaluated  by  quadrature. 

This  concludes  the  analytic  phase  of  the  ring-symmetric  CRW  approximation.  We  now  turn  to  the 
task  of  deriving  closed-form  expressions  approximating  the  flow  in  a  ring- symmetric  CRW  using  the 
corner  gradients  (v±0)a  . 


3.2      The   Matched   Extrapolation   Scheme 

We  seek  an  approximation  to  a  ring- symmetric  CRW  that  would  be  valid  up  to  several  corner  radii 
from  the  corner.  A  direct  substitution  of  the  C  *  corner  derivatives  (3.1-3)  in  a  Taylor  series 
truncated  after  the  first-order  term,  could  be  made.  The  result  would  be  an  expression  where  the 
local  deviation  of  the  ring-symmetric  CRW  from  the  corresponding  associate  CRW  is  proportional  to 
Ay  =  (y~~y  )  ,    i.e.,  a  'linear  extrapolation". 

It  turns  out  that  in  this  case  it  is  possible  to  formulate  a  "matched  extrapolation"  scheme. 
where  the  the  aforementioned  deviation  depends  nonlinearly  on  Ay.  which  would  be  consistent  to 
first-order  in  (Ay/y  )  with  the  corresponding  linear  extrapolation.  This  extrapolation  is  formulated 
to  match  exactly  the  flow  along  the  leading  characteristic  C  '  (pj).  The  advantage  of  this  improved 
scheme  is  that  it  maintains  better  accuracy  at  (Ay/y.)  =  0(1)  than  the  linear  extrapolation,  even 
though  either  scheme  is  formally  first-order  accurate  near  the  corner. 

The  implementation  of  the  matched  scheme  is  effected  through  a  change  in  flow  variables  from 
(v,0)  to  (f,0)  ,  where  f  is  the  streamtube  area  ratio  function  [8]  ,  linked  to  the  Prandtl-Meyer 
function  and  to  the  flow  Mach  number  through  the  isentropic  relations  : 
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-1    f  r  2-,i  (y+  l)/2(y-l) 

f  =  M     {(2/(Y+l))[l  +  ((Y-l)/2)M  ]) 

df/f  =  (M2-l)l/2dv  (3.2-1) 

-1   r  2rl   t       2  xl/2 

dv  =  M     [l  +  ((Y-l)/2)M  ]     (M-l)      dM 


We  note  that  relations  (3.2-1)  hold  in  any  steady  flow  where  entropy  is  uniformly  constant 
(homentropic  flow).  Thus,  the  thermodynamic  flow  variables  M  ,  f  and  v  can  be  used 
interchangeably.  The  replacement  of  v  by  f  is  motivated  by  reasons  which  will  be  made  clear  in 
the  sequel. 

Let  us  find  the  expression  for  fa(0,|3)  .  This  is  done  by  first  eliminating  9a(0.P)  by  summing 
(v  +  G)a  and  (v-0)a  in  (3.1-3),  obtaining  an  expression  for  va(0,P)  .  Then,  using  the  isentropic 
relations  (3.2-1)  to  replace   va  by  fa  and  using  (3.1-5)  for  ila(0.P)   and  <;an(0,p)   we  finally  get  : 


fu(0,P)  /  f(0,P)  =  6(0.p)  ya(0.P)  /  y(0,P) 

5(0,p)  =  (1/2)  {(p2-  i)I/2/  JJ,  na(0,P)  sinv|/(0,P)}    + 

(1/2)  {(P  -  l)    "  sin6(0.p)  /  p  sinM/<0\P)}    +    ((y+  l)/2(3-y))  H(P)  (3.2-2) 

H(P)  -  ((3-Y)/4)  {(P2-1)3/4  [l  +  ((y-l)/2)p2]'my~l)  I  siniK(O.p)}  * 

,•     r  ?i(2-Y)/2(7-l)  /     ,         x-5/4      , 

j     [l  +  ((Y-l)/2)m23  (m--l)       m2  sin6(0.m)  dm 

P, 

where  the  specially  defined  function  H(p)  has  been  normalized  by  a  preceding  factor  so  that 
H(oo)  =  1  .  A  typical  case  of  H(p)  is  shown  in  Figure  3-1.  By  pre-computing  (and  storing)  values 
of  H(p)  as  function  of  1/p  at  a  sequence  of  points  (typically  50  points),  H(p)  can  readily  be 
evaluated  by  linear  interpolation. 
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Consider  equations  (3.2-2)  at  the  leading  characteristic  C     (P<)  .     It  is  readily  verified  that 
H(Pj)  =  0  and  6(0,0,)  =  1  .   Consequently,  equation  (3.2-2)  assumes  the  following  form  at  the  leading 

i' 


characteristic  C     (P.)  : 


fa(0,P1)/f(0,p1)  =  ya(0,p1)/y(0,P1)  (3.2-3) 

Now,  along    C     (P,)    the  relation  between    f     and    ya    can  be  derived  directly  from  the  the 

equation  along   C      for  (v-8)a   (equation  (3.1-1)),  since  along   C     (P,)    0a  =  O,   and   v     can  be 

replaced   by   an   "isentropically   equivalent"   expression   containing     f    and     f      according   to    the 

isentropic  relation  (3.2-1).    The  result  is  the  following  differential  relation  which  holds  along  the 

entire    leading  characteristic  C     (p.)   of  the  ring- symmetric  CRW  : 

fa(o,p1)/«[oip1)  =  ya(o,p,)/y(atpl)  (3.2-4) 

This  differential  relation  can  be  readily  integrated,  yielding  f(a,P,)/y(a.p,)  =  constant  .  Indeed. 
this  result  is  consistent  with  the  very  definition  of  f  as  an  area  ratio  function,  since  the  flow  crossing 
the  leading  characteristic  C     (Pj)   diverges  with  cylindrical  symmetry. 

The  analogy  between  (3.2-3)  and  (3.2-4)  is  appealing.  Can  the  first  equation  of  (3.2-2)  be  recast  in 
a  form  that  would  yield  the  exact  relation  (3.2-4)  along  C  '  (P,)  .  rather  than  the  "linearized" 
relation  (3.2-3)  ?  Let  us  formally  rewrite  the  first  equation  of  (3.2-2)  as  if  it  were  valid  for  any 
a  £  0  ;   not  just  for   a  =  0  : 

fa(a,P)/f(a,P)  =  5(a,P)  [ya(a,P)/y(a,P)|  (3.2-5) 

Since  6(«.P)  is  not  known  to  us  (indeed,  it  is  defined  by  the  preceding  equation),  we  introduce 
the  following  assumption  which  constitutes  the  key  approximation  of  our  matched  extrapolation 
scheme  : 

8(a,P)  -  5(0.p)  (3.2-6) 

Using   this   approximation,   equation  (3.2-5)  is   readily   integrated,   yielding  a   relation   between 
f  and   y   that  does  not  explicitly  involve  the  function  y(a.p)  : 

f(a,p)  =  f(0,P)  [y(a,P)/y(0,P)]5(0,Pj  (3.2-7) 
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This  is  the  matched  extrapolation  of  the  area  function  f  ;  it  is  consistent  with  the  linear 
extrapolation  of  f  near  the  corner  as  would  be  derived  from  (3.2-2),  and  it  is  exact  along  the  leading 
characteristic  C  (Pj)  .  Does  this  matching  assure  a  reasonable  level  of  accuracy  throughout  the 
ring-symmetric  CRW  ?  We  defer  the  consideration  of  this  question  to  chapter  4,  where  the  subject  of 
accuracy  will  be  discussed  in  conjunction  with  comparison  to  finite-difference  computations  using  the 
code  JET. 

Let  us  examine  the  range  of  variation  of  6(0. P)  (Figure  3-2).  Along  the  leading  characteristic 
C  (pj)  we  have  5(0, pj)  =  1  ,  which  corresponds  to  a  cylindrically  divergent  flow.  At  high  values 
of  p  the  power  6(0, p)  increases  almost  monotonically  to  6(0, oo)  =  2/(3-y)  .  (It  actually  attains  a 
maximum  value  slightly  higher  than  6(0,  so)  --  see  Figure  3-2).  Assuming  that  a  dilute  gas  has  a 
specific  heats  ratio  y  ^  (5/3)  ,  the  maximum  value  of  6(0,^0)  for  any  dilute  gas  is 
6(0,^0)  ^  (3/2)  ,  which  corresponds  to  a  power  no  larger  than  midway  between  the  cylindrical 
power  (6=  1)  and  the  spherical  power  (6  =  2). 

The  power-law  form  of  approximation  (3.2-7)  to  streamtube  area  ratio  in  a  ring-symmetric  CRW 
suggests  the  following  geometrical  interpretation. 

(a)  Consider  the  variation  of  area  ratio  along  a  particular  streamtube.  It  consists  of  the  product  of 
two  factors  :  the  'planar  CRW"  area  ratio  f(0.P)  and  a  radial  divergence  factor,  where  the 
radial  power  is  function  only  of  P  .  There  is  some  form  of  separation  of  variables  in  (3.2-7). 
1(0. p)  and  0(0. p)  depend  only  on  P  which  designates  a  characteristic  line  C  '  (p)  .  The 
radial  divergence  factor  involves  y  which  designates  points  on  C"r(p).  This  is  a  direct 
consequence  of  our  key  assumption  (3.2-6). 

(b)  Consider  the  spacing  between  streamlines  in  the  (x.y)  plane  of  a  ring- symmetric  CRW.  They 
are  stretched  by  a  factor  (y/y  )  ^°'^1  relative  to  corresponding  streamlines  in  an  otherwise 
identical  planar  CRW.  The  stretching  power  vanishes  at  the  leading  characteristic  C'Np,)  and 
increases  to  an  asymptotic  value  of  (y-l)/(3-Y)   as   p   increases  to  infinity. 

(c)  At  very7  large  distance  from  the  corner,  one  would  expect  streamtube  cross-section  area  to 
exhibit  spherical  divergence,  i.e,  for  a  certain  streamtube  f  should  increase  as  y2  .  regardless 
of  the  value  of  y  .  The  fact  that  our  power-law  approximation  does  not  exhibit  these  features 
indicates  that  it  is  an  intermediate-range  fit  to  the  radially  divergent  CRW  flow.  Can  a 
matching  scheme  between  the  intermediate-range  and  the  far-range  flow  regimes  be  formulated? 
This  question  is  presently  an  open  one. 
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Any  flow  variable  linked  isentropically  to  f(a,p)  can  now  be  evaluated  from  the  approximation 
(3.2-7).  The  "kinematic"  variable  6(a,P)  ,  however,  is  not  related  to  f(a.P)  in  this  manner,  so  its 
approximation  is  still  pending. 

A  matched  extrapolation  scheme  for  G(a,p)  is  now  formulated  in  a  manner  consistent  with  the 
scheme  for  f(a,P)  .  This  is  done  by  regarding  equation  (3.1-1)  along  C  (P)  as  defining  a 
differential  relation  between  G(a.P)  and  y(a.P)  ,  since  the  differential  relation  between  v(a,P)  and 
y(a,P)  can  be  eliminated  by  virtue  of  (3.2-5)  and  the  isentropic  relation  between  v  and  f  given  in 
(3.2-1).  The  result  is  a  differential  relation  for  G(a,P)  along  C  (P)  analogous  to  the  relation 
(3.2-5)  for  f(a,p): 

9a(a,p)  =  c(«,P)  [ya(a,p)/y(a,p)] 

(3.2-8) 
c(a,P)=  6(a,p)  tann(a,p)  -  sinn(a,P)  sin0(a,p)  /  sin\}/(a,P) 

Now  we  invoke  the  same  argument  that  lead  to  the  key  approximation  (3.2-6),  namely,  that  £(a.P) 
may  be  replaced  by  its  associate  CRW  value  £(0,P)  .  Hence,  the  matched  extrapolation  for  G(a.p) 
is  given  by  : 

0(a.P)  =  9(0.p)  +  £(0.P)  ln[y(a.p)/y(0.p)j 

(3.2-9) 
£(0.p)  =  6(0. p)  tann(O.p)  -  sin^O.p)  sinG(O.p)  /  siny(O.P) 

Let  us  examine  the  function  £(0.p)  (Figure  3-3).  We  observe  that  £(0,p,)  =  0  in  accordance  with 
the  boundary  condition  at  the  leading  characteristic  C  (Pj)  :  8)X(a.P)=0.  (Since  by  (3.2-8) 
£(a,Pj)  =  0  for  any  a  ,  this  boundary  condition  cannot  be  considered  as  a  guideline  in  improving  the 
6(a.P)  extrapolation  as  was  the  case  with  the  f(a.p)  extrapolation).  At  large  p:  5(0.p)->2/(3-y)  , 
so  that  p£(0.P)— »("/-  l)/(3-y)  .  Consequently,  the  effect  of  ring  symmetry  is  to  deflect  streamlines 
toward  the  radial  direction  (y)  relative  to  their  associate  CRW  direction,  by  an  amount  that  decreases 
in  the  hypersonic  portion  of  the  CRW.  as  the  inverse  of  the  associate  Mach  number. 

An  unexpected  difficulty  arises  when  the  logarithmic  extrapolation  (3.2-9)  for  6(a.P)  is  considered 
for  large  y(a,P)  .  At  some  sufficiently  large  v(a.p)  ,  G(a.P)  will  exceed  nil  ,  which  seems 
physically  unrealistic.  A  possible  remedy  would  be  to  replace  (3.2-9)  by  a  modified  extrapolation 
obtained    as    follows.     Multiply   (3.2-8)   by     tanG(a.P)  ,     and   invoke    the    modified   approximating 
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assumption       that        c(a,P)  tan0(a,p)        may       be       replaced       by        £(0,P)  tanG(0,p)  .        Since 
d(tanz)  =  -  d(ln(cosz))  ,  the  result  is  the  following  power-law  extrapolation  : 

r  i-K(0,p) 

cos6(a,p)  =  cosG(0,P)  [y(a,P)/y(0,P)J 

(3.2-10) 
K(0,P)  =  £(0?P)  tan9(0,P) 

We  notice  that  by  this  power-law  approximation  6(a,p)  cannot  exceed  n/2  .  Also,  in  applying 
(3.2-10)  care  should  be  exercised  near  p  =  P1?  since  K(0,pj)  is  singular.  However,  a  simple  analysis 
shows  that  K(0,p)  approaches  a  finite  limit  at  pL  ,  which  can  be  expressed  in  closed  form,  and  since 
cos0(O,Pj)  =  O  the  boundary  condition  cosG(ct,pj)  =  0  is  fulfilled  regardless  of  the  value  of  K(0,Pj)  . 
The  power  K(0,P)   for  a  typical  case  is  shown  in  Figure  3-4. 

One  more  variable  is  needed  to  complete  the  approximation  to  the  ring- symmetric  CRW  :  the 
C     (p)  characteristic  direction  \|/(ot,P)  •   It  is  simply  given  by  : 

V/(a,P)  =  G(a,P)  +  H(a,p)  (3.2-11) 

where  0(a.p)  is  given  by  either  (3.2-9)  or  (3.2-10).  The  Mach  angle  fi(a.p)  is  obtained  from 
f(a,P)   through  the  isentropic  relation  (3.2-1),  where   f(a.P)   is  given  by  the  approximation  (3.2-7). 

This  concludes  the  specification  of  a  ring- symmetric  CRW  through  matched  extrapolation  of  the 
associate  CRW  .  To  obtain  the  flow  field  of  a  particular  ring-symmetric  CRW,  one  integrates  the 
relation  defining  C  (P)  ,  i.e..  dx  =  cot\|/dy,  while  evaluating  the  flow  variables  using  the 
aforementioned  approximations.  The  result  is  flow  field  information  at  a  series  of  points  (x.y)  along 
C     (P)   characteristic  lines. 

There  remains  one  more  approximation  task  :  the  "inverse  problem",  which  is  stated  as  follows. 
Given  a  point  (x«,y0)  in  the  ring- symmetric  CRW,  trace  'inversely"  the  characteristic  line  through  it, 
and  hence  find  the  flow  variables  solely  at  (X0,y0)  ,  using  the  matched  extrapolation  scheme.  An 
efficient  method  for  resolving  the  inverse  problem  is  described  in  section  3.3  below. 
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3.3      The   Inverse   Problem 

Consider  the  following  application  of  the  present  approximation  to  a  ring- symmetric  CRW.  A 
physical  model  is  related  to  the  flow  in  a  ring- symmetric  CRW.  The  computation  scheme  of  that 
model  requires  flow  variables  at  a  large  number  of  points  (x0,y0)  within  the  flow  field.  The 
implementation  of  such  a  model  calls  for  an  efficient  algorithm  capable  of  providing  flow  variables  at 
points  (xQ.y0)  ,  without  resorting  to  computation  of  the  flow  at  a  surrounding  cluster  of  points,  or 
along  a  set  of  adjacent  characteristic  lines. 

It  is  observed  that  through  all  but  the  low  supersonic  portion  of  a  ring-symmetric  CRW,  the   C 
characteristics  are  nearly  straight  lines.    This  suggests  the  following  approximation  procedure  for  the 
inverse  problem.   Assume  as  a  first  guess  that   C     (P)   through  point   (xQ,y0)   is  a  straight  line;   then 
correct  the  resulting  value  of  P    by  finding  the  small  deviation  of  C     (p)   from  a  straight  line,  using 
a  linearized  approximation  to  vj/(a,P)  . 

To  do  that  we  need  a  closed-form  expression  for  V|/(a,P)  -  vj/(0,P)  .  The  implicit  definition  (3.2-1 1) 
will  not  be  adequate.  Instead,  we  seek,  an  extrapolation  scheme  formulated  directly  for  y(a.p)  . 
Following  a  procedure  analogous  to  that  by  which  the  f(a.P)  and  6(a.P)  approximations  were 
obtained,  we  get  : 

i|/a(a,p)  =  X.(a,P)  [ya(a,P)/y(a,p>] 

A(a.P)  =  e(a.p)  -  6<a.P)  [l  +  ((Y-l)/2)M"]  (m~-i)~  (3.3-1) 

M  =  M(a.p) 


Again  we  invoke  the  argument  that  (3.3-1)  is  exact,  but  it  involves  the  unknown  functions  G(a.B)  . 
6(a,P)  and  M(a.p)  .  so  that  they  ought  to  be  approximated  by  their  respective  associate  functions. 
i.e.,  S(0,p)  .  6(0. p)  and  M(0.p)  .  (  We  might  have  retained  the  Mach  number  as  M(a.P)  and 
approximate  it  through  (3.2-7);  however,  that  would  have  made  it  impossible  to  obtain  a  closed- 
form  expression  for  the  ensuing   y  -  integration).   Thus,  replacing  X(a.P)   by  X(O.P)  ,   we  get  : 


19 


V(a,P)  -  v(0,P)  =  X(0,P)  ln[y(a,P)/y(0,p)] 

X(0,p)  =  c(0,p)  -  8(0,P)  [1  +  ((Y-1)/2)P2]  (p2"l)"3/2 


The  unknown  Pn  is  now  to  be  determined  by  solving  the  equation  : 


vo 


(3.3-2) 


x0  =  xc  +      J      coty(a,p0)  dy(a,p0)  (3.3-3) 

jc 


To  get  an  explicit  solution  for  p0  from  (3.3-3),  we  "linearize"  this  equation  by  letting  : 
coti^p,,)  *  coti|/(0,P0)  -  [  \|/(a,P)  -  v(0,P0)  ]  /  sin2V|/(0,P0)  (3.3-4) 

We  substitute  this  approximation  in  (3.3-3)  and  we  use  the  expression  (3.3-2)  for  y(a.P)  .   The   y 
integration  can  then  be  performed  in  closed  form,  yielding  the  following  equation  for  p0  : 

cot\|/(0.p0>  =  (x0-xc)/(y0-yc)  + 

I  y0  ln(y0/yc)  /  (y0~yc)  "  [  1  [Mo.P<p  /  sin2v(o,p0)]  (3.3-5) 


When  X(0.po)  is  sufficiently  small,  (3.3-5)  can  be  solved  for  P0  by  repeatedly  computing  p()  on 
the  left  side  using  a  former  iteration  value  of  p0  on  the  right  side.  In  this  procedure,  we  also  use  the 
associate  CRW  relation  between  vy(0,p0)  and  P0  (see  Appendix  A).  For  the  initial  iteration  we  set 
X(0,p0)  =  0. 

The  accuracy  o[  solving  the  inverse  problem  is  not  as  good  as  that  of  the  direct  matched 
extrapolation.  Additional  errors  are  generated  by  the  linearization  (3.3-4)  and  also  by  the 
approximation  (3.3-2)  to  Vj/(a,p)  which  is  inferior  to  (3.2-11).  Thus,  along  the  leading  characteristic 
C  '  (pj)  the  direct  approximation  becomes  exact,  but  the  inverse  procedure  remains  an 
approximation. 
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The  coefficient  X(0,p)  is  usually  fairly  small  (Figure  3-5).  For  large  P  the  following  asymptotic 
expression  is  readily  found  from  (3.3-2)  :  pX(0,P)-»0  ,  so  that  the  accuracy  of  the  inverse  procedure 
improves  as  (x^y^  moves  into  the  hypersonic  portion  of  the  ring-symmetric  CRW.  If  this  accuracy 
is  inadequate  when  (xQ,y0)  is  near  C  (P,)  ,  we  can  revert  to  solving  (3.3-3)  for  P0  by  some 
iterative  scheme,  performing  the  integration  along  C  (P0)  numerically,  while  using  the  better 
approximation  (3.2-11)  for  y(a,P0)  .  Naturally,  this  procedure  would  entail  a  much  higher 
expenditure  of  CPU  time  than  the  approximation  (3.3-5). 

For  the  sake  of  completeness,  there  is  one  more  point  to  discuss  in  regards  to  the  inverse  problem. 
When  the  iterative  procedure  of  Eq.  (3.3-5)  involves  values  of  P0  sufficiently  close  to  the  value  Pj 
of  the  leading  characteristic,  a  difficulty  may  arise.  The  reason  is  that  while  the  leading  characteristic 
is  exactly  replicated  by  the  matched  approximation,  the  "linearized"  approximation  to  the  leading 
characteristic  (obtained  by  computing  XQ  as  function  of  y0  from  Eq.  (3.3-5)  with  P0=P|  ),  curves 
towards  the  interior  of  the  CRW.  This  situation  is  depicted  in  Figure  3-6.  Thus,  the  iteration 
procedure  of  Eq.  (3.3-5)  would  come  up  with  a  tentative  solution  having  P0  <  Pj  for  points  (xQ,y0) 
located  between  the  exact  and  approximate  leading  characteristics  (shaded  area  in  Figure  3-6),  which 
is  inconsistent  since  the  range  of  acceptable  values  within  the  CRW  is  P0  >  Pr  The  remedy  is 
simply  to  set  P0  =  P T  as  the  approximate  solution  to  any  point  (x0,y0)  within  that  "shaded  area". 
To  illustrate  the  magnitude  of  error  introduced  by  this  simplification,  we  show  in  Figure  3-6  the 
characteristic  corresponding  to  p  =  4.1  (versus  P,=4  ).  This  difficulty,  however,  is  probably  of  no 
concern  to  ambient  scattering  applications,  since  ambient  molecules  are  typically  stopped  at  the  outer 
fringes  (hypersonic  portion)  of  the  ring-symmetric  CRW. 
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Fisure   3-1. 


The  Integral  Function  H(P)   as  Function  of  Inverse  Mach  Number 
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Fieure   3-2. 


Radial  Power  for  Area  Ratio  as  Function  of  Inverse  Mach  Number 
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Figure   3-3.  Coefficient  of  Logarithmic  Variation  in  Flow  Angle  as  Function  of  Inverse  Mach 

Number 
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Figure  3-4. 


Radial  Power  for  Cosme  of  Flow  Ansle  as  Function  of  Inverse  Mach  Number 
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Figure   3-5.  Coefficient  of  Logarithmic  Variation  in  characteristic  Angle  as  Function  of  Inverse 

Mach  Number 
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between  Exact  and  Linearized  Leading  Characteristics 
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4.        RESULTS   AND   DISCUSSION 

We  now  present  computational  results  intended  to  exhibit  some  relevant  features  of  ring- symmetric 
supersonic  jets,  and  to  demonstrate  the  nature  of  our  matched  approximation  to  the  CRW  region  of 
the  ring-jet  flow  field.  A  typical  case  of  HF/DF  laser  exhaust  was  chosen  [1,10]  for  these  sample 
computations.   It  is  specified  in  Table  4-1  below. 


Table  4-1.  Typical  Operating  Conditions  of  HF/DF  Laser  Exhaust 


Mole  fractions 


[H  ]    =  .091     [HF]  =  .091     [H2]  =  .104 
[DF]  =  .135    [He]  =  .579 


Average  molecular  weight 


7.27 


Specific  heats  ratio 


1.54 


Stagnation  temnerature  and  densitv  2300    (K)    and    .0075    fks/nr) 


Exit  Mach  number 


4.0 


Molecular  diameter  (hard  spheres) 


2.5xl0"10    (m) 


Spacecraft  diameter 


Aperture  of  rinsz-nozzle 


2.5    (m) 
1.0    (m) 
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4.1      Finite   Difference   Computation  of  Ring-Symmetric   CRW 

An  accurate  finite  difference  computation  of  the  entire  ring- symmetric  CRW  flow  field  was 
performed  by  the  code  JET  utilizing  the  SI  MA  scheme  (Ch.  2).  A  brief  description  of  the 
computational  procedure  is  now  presented.  Referring  to  Figure  4-1,  the  y  marching  step  started  at 
0.01  (m)  and  was  subsequently  limited  to  half  the  step  determined  by  the  closest  forward  intersection 
of  pairs  of  C  characteristics  extended  from  adjacent  grid  points.  At  the  outset,  the  number  of 
lateral  (x)  grid  intervals  was  100,  of  which  60  were  in  the  CRW  and  40  spanned  the  nozzle  exit. 
SI  VIA  scheme  was  applied  to  all  CRW  grid  points  which  were  determined  by  a  fan  of  C 
characteristic  lines  indexed  from  k=  1  (leading)  to  k=61  (final).  These  fan  lines  were  chosen  to  have 
equally  spaced  values  of  the  Riemann  invariant  (v  -  0)  at  the  corner.  The  flow  at  all  remaining  grid 
points  was  computed  via  the  inverse  marching  scheme  (Ch.  2). 

As  the  marching  progressed,  the  leading  characteristic  (k=  1)  approached  the  mid-plane  of 
symmetry  (x=0)  and  it  became  necessary  to  switch  the  first  SIMA-computed  characteristic  to  k=  2. 
then  subsequently  to  k=3  and  so  on,  until  at  y  =  24  (m)  it  reached  the  index  k=19  (Figure  4-1). 
The  result  is  that  in  addition  to  the  region  designated  as  CRW  (Figure  2-2),  also  the  region  roughly 
coinciding  with  the  open  triangle  between  k=  19  on  the  left  and  the  reflected  characteristic  k=  I  on 
the  neht,  was  comDuted  via  the  SIM  A  scheme. 


4.2      Continuum   Breakdown 

In  a  source-like  flow  of  dense  gas  exhausting  into  vacuum,  a  breakdown  of  continuum  flow 
inevitably  takes  place  at  some  point  along  each  streamtube  [2].  In  a  planar  CRW  i  Prandtl-Mever 
flow),  the  breakdown  surface  approaches  asymptotically  a  particular  streamline  as  the  distance  from 
the  corner  increases  [3.10].  However,  in  a  ring-symmetric  CRW  the  breakdown  surface  curves 
towards  the  mid-plane  ot"  symmetry,  reaching  it  at  some  finite  radius.  Two  breakdown  surfaces  were 
traced  in  the  JET  computation,  corresponding  to  the  plausible  values  of  the  breakdown  parameter 
B=0.05  and  B=0.08  ;  they  are  shown  in  Figure  4-2.  It  is  evident  from  this  figure  that  most  of  the 
ring- symmetric  CRW  lies    outside     the  breakdown  surface. 

Moreover,  consider  ambient  scattering  by  molecules  entering  the  CRW  in  the  x  direction.  Most 
of  these  molecules  will  be  stopped  within  a  range  for  which  the  expected  number  of  collisions  with 
exhaust  molecules  is   T~  1  .    In  Figure  4-3  we  show  the  line    x  =  1    plotted  from  the  JET  computation 
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of  the  typical  laser  exhaust  (Table  4-1).  It  is  clear  that  virtually  all  ambient  scattering  will  take  place 
outside  the  breakdown  surface.  Our  model  for  estimating  the  molecular  backflow  induced  by  ambient 
scattering,  relies  on  the  presumption  that  continuum  flow  holds  throughout  the  rarefaction  fan,  even 
at  regions  of  near- vacuum.   How  can  that  approach  be  justified  in  view  of  the  preceding  discussion? 

Consider  a  spherical  source  flow  into  vacuum  [2]  ,  where  the  breakdown  surface  is  in  a  region  of 
hypersonic  flow.  A  special  situation  in  regards  to  continuum  breakdown  exists  in  this  idealized  case. 
The  breakdown  in  continuum  flow  is  first  manifested  by  the  failing  of  local  temperature  to  adhere  to 
the  isentropic  relation  with  local  density.  In  fact,  "temperature"  ceases  to  exist,  since  the  random 
component  of  molecular  velocity  does  no  longer  have  the  (isotropic)  Boltzmann  distribution.  By 
contrast,  conservation  of  mass  holds  regardless  of  continuum  flow  breakdown,  so  that  by  virtue  of  the 
spherical  symmetry,  density  decreases  as    1/R"  to  arbitrarily  large  R  . 

As  demonstrated  by  Bird's  Monte  Carlo  simulation  of  a  Prandtl-Meyer  flow  [3]  ,  an  analogous 
situation  exists  also  in  this  self-similar  flow  field  :  the  density  follows  quite  closely  the  continuum 
(self  similar)  solution,  even  in  regions  well  downstream  from  points  where  the  isentropic  relation 
between  temperature  and  density  has  broken  down.  Since  only  density  and  flow  velocity  are  needed 
in  order  to  estimate  ambient  scattering  (see  [1]  .  in  particular  the  "cold"  assumption),  the  foregoing 
discussion  provides  a  justification  for  the  use  of  continuum  flow  approximation  to  this  end. 


4.3      The   Matched   Approximation 

Consider  a  particular  characteristic  line  (k  =  47)  in  the  sample  case  computation  (Figure  4-1).  We 
focus  on  the  variation  of  the  three  thermodynamic  variables  f  (area  ratio).  M  I  Mach  number)  and  v 
(Prandtl-Meyer  function),  along  this  characteristic  line.  These  variables  were  computed  via  both  the 
SI  MA  scheme  (Ch.  2)  and  the  matched  approximation,  and  the  results  are  shown  in  Figures  4-4,  4-5 
and  4-6.  As  expected,  Mach  number  and  area  ratio  both  increase  along  the  characteristic  line,  due  to 
the  radial  divergence  of  the  flow.  How  well  is  this  variation  predicted  by  our  matched 
approximation? 

The  matched  approximation  for  thermodynamic  variables  is  the  radial  power-law  expression  for 
the  area  ratio  (3.2-7).  M  and  v  are  computed  from  f  via  the  standard  isentropic  relations.  The 
following  features  are  observed  regarding  the  nature  of  the  matched  approximation  in  this  sample 
case  (Figures  4-4,  4-5  and  4-6). 


30 


(a)  Both  SIMA  and  matched  approximation  curves  start  out  with  equal  corner  gradients.  This 
constitutes  a  mutual  validation  for  both  the  numerical  integration  scheme  and  the  analytic 
approximation,  in  the  vicinity  of  the  corner. 

(b)  The  approximated  f  (Figure  4-4)  is  undervalued,  resulting  in  a  correspondingly  undervalued 
Mach  number  (Figure  4-5).  Had  we  opted  for  a  linear  extrapolation  of  f  from  the  corner 
along  fan  characteristics,  it  is  evident  from  Figure  4-4  that  f  would  have  been  even  more 
undervalued,  and  so  would  be  the  Mach  number.  Matched  approximation  is  thus 
demonstrated  to  be  clearly  superior  in  accuracy  to  linear  extrapolation. 

(c)  What  would  be  the  quality  of  an  approximation  based  on  a  linear  extrapolation  of  v  ?  Using 
this  variable  would  be  consistent  with  the  fact  that  the  analytic  expression  of  fan-wise 
gradients  (3.1-3)  is  for  gradients  of  the  Riemann  invariants. 

From  Figure  4-6  it  is  evident  that  a  linearly  extrapolated  v  would  be  grossly  overvalued, 
while  the  matched  approximation  for  v  is  moderately  undervalued.  A  linear  extrapolation  of 
v  would  thus  result  in  a  grossly  overvalued  Mach  number.  (In  fact,  from  Figure  4-6  it  can  be 
estimated  that  this  would  result  in  M  reaching  infinity  around  y/y  =4  ).  This  observation 
is  not  inconsistent  with  (b)  above,  as  may  seem  at  first.  The  isentropic  relations  between  f  . 
M  and  v  are  highly  nonlinear,  so  that  the  quality  of  an  approximation  scheme  depends  on 
the  particular  thermodynamic  variable  chosen  to  be  approximated  (while  the  others  are 
evaluated  from  the  isentropic  relations).  The  observation  that  a  linear  extrapolation  of  v 
would  be  a  poor  cnoice  of  approximation  scheme,  thus  provides  yet  another  support  to  the 
choice  of  f  as  the  thermodynamic  variable  to  be  approximated. 

Can  any  statement  be  made  at  this  point  about  the  level  of  accuracy  of  matched  approximations 
to  ring-symmetric  CRW  ?  Discussing  this  question  is  contingent  upon  specifying  some  physical 
feature  of  the  flow  field  for  which  a  comparison  will  be  made.  Merely  computing  a  How  variable, 
such  as  Mach  number,  does  not  constitute  a  meaningful  test. 

For  the  application  which  motivated  the  present  study,  namely  ambient  scattering,  a  meaningful 
feature  of  the  How  field  would  be  the  cumulative  number  of  collisions  that  a  fast  moving  molecule 
entering  the  CRW  from  a  bordering  cavitation  region,  can  expect.  This  quantity  is  path-dependent, 
so  as  a  simplification  we  consider  only  rays  parallel  to  the  x  axis  (constant  y  ).  The  number  of 
collisions  expected  by  a  molecule  reaching  point   (x0,y0)   is  : 

t(x0,y0)  =      J      <jn(x,y0)dx  (4.3-1) 

x, 
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Where  n  is  the  number  density,  <r  is  the  collision  cross-section  and  xL  corresponds  to  the 
limiting  characteristic.  We  generally  assume  hard-sphere  collisions  so  that  <J  is  a  constant,  and  the 
"molecular  opacity"   t  is  essentially  a  mass  integral  multiplied  by  a. 

The  molecular  opacity  T(x0,y0)  was  computed  along  with  the  computation  of  the  flow  field  by 
the  code  JET.  Some  lines  of  constant  opacity  were  traced;  two  of  them  are  shown  in  Figure  4-3.  The 
opacity  in  the  CRW  region  was  also  computed  from  the  matched  approximation;  the  results  compare 
well  with  those  of  SI  MA  computations,  as  can  be  seen  in  Figures  4-7  and  4-8.  The  method  by  which 
t(x0.y0)  was  evaluated  from  the  matched  approximation  is  not  a  numerical  integration.  Rather,  it  is 
a  closed-form  approximation  which  we  now  proceed  to  describe. 

Referring  to  Figure  4-9,  the  law  of  conservation  of  mass  along  a  particular  streamtube  is  used  to 
establish  the  following  relation  between  flow  variables  at  points  S  on  the  characteristic  line  passing 
through  (x0,y0)   and  points  X   on  the  segment   y=y0: 

n(S)  u(S)  sinfl(S)  y(S)  AS  =  n(x)  u(x)  sinG(x)  y0  Ax  (4.3-2) 

The  opacity  integration  (4.3-1)  can  thus  be  replaced  by  an  integration  along  the  characteristic  line 
from  (xc,yc)  to  (X^y^)  .  Now,  some  further  simplifications  are  introduced,  enabling  a  closed-form 
integration  for  the  opacity.  First  we  assume  that  the  entire  region  of  interest  is  one  of  verv  msih 
Mach  number,  so  that  I  +  ((y-  1)/2)M  «  ((y-l)/2)M"  .  Powers  of  M  can  be  approximately 
replaced  by  appropriate  powers  of  f  .  We  also  assume  that  \\)  and  0  do  not  vary  greatly  through 
the  region  of  interest,  so  that  both  can  be  reasonably  well  approximated  by  their  corner  values. 
Using  the  power  law  (3.2-7),  the  opacity  integration  is  then  readily  performed,  yielding  the  following 
expression  : 

t(x0,y0)  =  G<0.p())(y0/y/  [(y0/yc)2^+1^2  -  l]  /  [2 - 6(7  +  l)/2] 

G<0,P0)  =  (7  n(0.(30)  yc  /  M(0.p0)  sin0(O.po)  sinn/iO.fy  (4.3-3) 

6  =  6(0.p0)  (see  Eq.(3.2-2)  above) 

The  opacity  distributions  shown  in  Figures  4-7  and  4-8  were  computed  from  this  expression.  They 
compare  favorably  with  SI. MA  opacities,  demonstrating  that  for  our  purpose  the  matched 
approximation  is  reasonably  accurate.     In  these  JET  computations  the  opacity  at  points   on  the 
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boundary  (final)  characteristic  line  was  estimated  from  Eq.  (4.3-3).  The  value  of  the  boundary' 
opacity  represnts  the  effective  thickness  of  the  "missing"  layer  of  fluid  due  to  terminating  the 
numerical  computation  at  the  characteristic  having  M(0,p)  =  Mf .  The  maximum  value  of  this 
boundary  opacity  in  the  typical  case  was  found  to  be  about  0.16  (at  y/yc  =  3.2  as  in  Fig.  4-7).  This 
relatively  low  value  indicates  that  as  far  as  the  interaction  with  invading  ambient  molecules  is 
concerned,  the  final  Mach  number  Mf=  34  is  a  reasonable  substitute  for  Mf  =  °o  . 

We  observe  that  Eq. (4.3-3)  predicts  that  opacity  generally  peaks  at  some  point  along  each 
characteristic  line  (due  to  its  particular  dependence  on  y  ),  which  in  the  case  shown  in  Figure  4-8,  is 
confirmed  by  the  SIMA  opacity  curve.  An  interesting  feature  of  the  expression  (4.3-3)  is  that  it  is 
given  in  a  separation-of-variables  form,  the  variables  being  P0  and  y0  that  uniquely  designate  a 
point  in  the  CRW.  We  believe  that  this  approach  can  be  adapted  to  other  physical  problems  calling 
for  opacity  integrals  in  a  ring- symmetric  CRW  . 
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Fisure  4-1.  Characteristics  in  Rine- Symmetric  Jet 
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Figure  4-2.  Continuum  Breakdown  Surfaces  in  a  Ring-Symmetric  Jet 
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Fieure  4-3. 


Constant  Molecular  Opacity  (t)  Lines  in  a  Ring- Symmetric  Jet 
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Variation  of  Area  Ratio  alone  Characteristic  Line  k  =  47 
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Figure  4-5. 


Variation  of  Mach  Number  along  Characteristic  Line  k  =  47 
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Figure  4-6.  Variation  of  Prandtl-Meyer  Function  along  Characteristic  Line  k=47 
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Figure  4-7.  Variation  of  Molecular  Opacity   I   along  a  Line  Parallel  to   x  axis 
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Figure  4-8.  Variation  of  Molecular  Opacity   T   along  a  Characteristic  Line  (k=55) 
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Figure   4-9. 


Schematic  Description  of  Mass  Conservation  Law  Enabling  Transformation  of 
Opacity  Integral  to  Characteristic  Line  Passing  through  (X0,y0) 
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APPENDIX  A.    PRANDTL  -  MEYER  FLOW 

Denoting  by    P    the  Mach  number  of  the  fan  characteristic,  the  Prandtl- Meyer  function  [8]    is 
given  by  : 

1/2  r,    2        >J/2       l/2i  IV    2        \l/2l 

v(P)  =  T       arctan[(p   -  l)  '  /r     J  -  arctanl(p   -  l)     J 

(A-l) 

r  =  (y+i)/(y-d 


In  a  C  wave,  the  Riemann  invariant  (v  +  G)  is  uniformly  constant.  This  leads  (after  some 
manipulation)  to  the  following  relations  between  the  Mach  number  P  and  the  C  characteristic 
angle  \\i  : 


V(P)  =  VL  +  r1/2arctan[r1/2/(p2-l),/2] 

(A-2) 
p2  =  1  +  r/tan2[(v-i|/L)/rI/2] 


Where  the  characteristic  ande  and  Macn  ansie  are  aiven  bv 


V(P)  =  9(p)  +  wp) 

(A-3) 
2        xI/2. 


M(P)  =  arcsin(p)  =  rc/2  -  arctan[(p~-  l)    "j 


We  now  seek  an  expression  for  the  distance  rj  of  a  point  on  the  transverse  (C  )  characteristic 
from  the  corner.  The  following  differential  relation  is  derived  from  flow  field  geometry  (Figure  A-l) 
in  a  planar  CRW  : 


r\  d\j/  -   -  tan(2n)  drj  (A-4) 
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In  order  to  replace    Ay    by   Ap  ,   we  use  the  following  expression  derived  by  differentiating  the 
relation  vji(P)   given  in  Eq.(A-5)  with  respect  to   P  : 

Vr  -  -  <(Y  +  0/2)  P  (P2-  l)"1/2  [l  +  ((Y-l)/2)p2r'  (A-5) 


Using  (A-5),  we  eliminate  Av|/   from  (A-4),  getting  : 


(lnn)p  =  ((Y  +  0/4)  P  (P2  -  2)  (P2  -  I)"'  [  1  +  ((y-  1)/2)P2]"1  ( A-6) 


The  right  hand  side  of  (A-6)  is  decomposed  into  two  partial  fractions,  resulting  in  an  integral 
consisting  of  the  sum  of  two  log  functions.  Thus,  the  dependence  of  v\  on  P  is  the  following 
product  of  two  factors,  normalized  to  yield   t|(Pi)  at  tne  leading  characteristic  : 

TUP)  =  rtfPj)  {[l  +  ((Y-l)/2}p23/{l  +  <(7-I)/2)p12]}r//2(Y"'      {(p,:-l)/(p2-l)}1/4  (A-7) 


Since     r|(p)     is    identical    with    the    definition    (3.1-2)    of  the    characteristic    variable     a  .    by 
differentiating  (A-7)  with  respect  to   a  with  P  =  constant,  we  get  Eq.i  3.1-5 1  for  11    . 

Now  we  turn  to  the  other  geometrical  function  needed  in  the  analysis  of  the  ring-symmetric  CRW 
flow  at  the  corner.  The  following  geometrical  relationship  is  established  by  inspection  from  Figure 
A-l  : 

A^  =  Art  /  cos(2n)  (A-8) 

Where  the  t  and  rj  increments  correspond  to  an  increment  in  characteristic  angle  \\i  .  or  to  an 
implied  increment  in  Mach  number  P  .  When  Ar|  is  replaced  by  AP  using  (A-6)  and  when  also 
cos(2^i)  is  replaced  by  a  p  expression,  we  get  a  expression  for  £^  in  terms  of  rj  and  p  .  The  a 
derivative  of  this  last  expression  is  Eq.(3.1-5)  for   ^ao  . 
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Geometry  of  Characteristics  and  Streamlines  in  a  Prandtl-Mever  Flow 
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